Defect stability in phase-field crystal models: Stacking faults and partial dislocations 
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The primary factors controlling defect stability in phase- field crystal (PFC) models are examined, 
with illustrative examples involving several existing variations of the model. Guidelines are presented 
for constructing models with stable defect structures that maintain high numerical efficiency The 
general framework combines both long-range elastic fields and basic features of atomic-level core 
structures, with defect dynamics operable over diffusive time scales. Fundamental elements of the 
resulting defect physics are characterized for the case of fee crystals. Stacking faults and split 
Shockley partial dislocations are stabilized for the first time within the PFC formalism, and various 
properties of associated defect structures are characterized. These include the dissociation width 
of perfect edge and screw dislocations, the effect of applied stresses on dissociation, Peierls strains 
for glide, and dynamic contraction of gliding pairs of partials. Our results in general are shown to 
compare favorably with continuum elastic theories and experimental findings. 
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I. INTRODUCTION 

Structural kinetics in crystalline solids are driven het- 
erogeneously at the atomic level by localized defects, 
which in turn drive mesoscopic and macroscopic phenom- 
ena such as structural phase transformations, fracture, 
and other forms of plastic flow. A complete description of 
such processes therefore requires a multiscale approach. 
Existing modeling methods typically operate either ex- 
clusively on atomic scales or on meso- and macroscopic 
scales. Phase-field crystal models on the other hand pro- 
vide a framework that combines atomic length scales and 
mesoscopic / diffusive time scales^^, with the potential to 
reach mesoscopic lengths through systematic multi-scale 
expansion method*^. 

The PFC approach naturally incorporates elasticity, 
plasticity, and effects of local crystal orientation into a 
relatively simple atomic-level continuum theory^. The 
literature published to date demonstrates that such a 
formulation, with only a few minimal ingredients, gives 
rise to a broad range of physics associated with diffu- 
sive nonequilibrium processes in liquid and solid systems. 
One of its strengths in terms of describing periodic sys- 
tems is the wealth of inherent defect structures that au- 
tomatically emerge from the basic free energy functional. 
This permits the study of crystalline defects within a de- 
scription that captures both long-range elastic fields and 
basic features of atomic-level core structures. In addi- 
tion to fundamental, local defect properties, the role of 
these defects in dynamic materials phenomena that op- 
erate over long, diffusive time scales can be examined. 
Therefore, defect stability, the core structure of stable 
defect configurations, and the dynamics and interactions 
of various defect structures within the PFC description 



all become central issues as the method is advanced fur- 
ther into solid-state phenomena. However, in light of the 
approximations inherent to the approach, its limitations 
need to be adequately realized and understood as well. 

The basic properties of perfect PFC crystals are rela- 
tively well understood, but the various classes and types 
of defects relevant to each given lattice symmetry require 
further study if material-specific applications are to be 
pursued. Perfect dislocations and simple grain bound- 
aries in 2D triangular and 3D bec crystals have been ex- 
amined and are known to be topologically stable under 
typical conditions of smal l local strain and low thermo- 
dynamic driving force^S^HUl Certain other defect struc- 
tures, notably stacking faults in close-packed 3D crystals, 
have inherently lower topological stability, and in such 
cases the proper balance between crystal stability and 
defect stability in the model formulation becomes more 
restrictive and difficult to achieve. It is also not neces- 
sarily clear which of the now many versions of PFC are 
best suited for describing defect-mediated processes. 

The aims of the first part of this article are to char- 
acterize the nature of defect stability in PFC models in 
terms of a few general model features, and to examine 
how various versions compare in this respect. We will 
show that fine control of crystal structure often leads to 
reduced defect stability. A balance must therefore be 
found that sufficiently favors crystalline order yet does 
not destabilize relatively unprotected defect structures 
such as planar faults. A few ways to achieve this balance 
are identified and compared critically. The most efficient 
of these approaches is then, in the second part of this 
article, applied to the case of fee crystals, in which the 
stability of stacking faults plays a central role. A range 
of defect properties and behaviors consistent with theory 
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and/or experiment are shown to naturally emerge from 
what is still a very simple set of equations. Potential 
implications for atomistic studies of plastic deformation 
involving slow, diffusive processes are discussed, and ini- 
tial results concerning a few processes of interest are de- 
scribed. More detailed findings relative to these will be 
outlined further in a future publication. 

PFC models can be viewed as simplifie d versi ons of 
classical density functional theory (CDFT^SEB. The 
Ramakrishnan-Yussouff (RY) CDFT^ provides a useful 
reference point, with a free energy functional given by an 
expansion around the liquid state correlation functions, 



applicable to Eq. ([!]), is 
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With imposed density conservation this equation pro- 
vides an accelerated path to local energy minima 2 -^. All 
simulations in this study were performed in 3D using 
pseudo-spectral algorithms and periodic boundary con- 
ditions. Those that employed Eqs. ^ or Q used semi- 
implicit time stepping, while those that employed Eq. ([5| 
used explicit time stepping. 



df \p(r) In (p(r)/pe) - Sp(f}} 



df df 2 Sp{f)C 2 {f,f 2 )Sp{f 2 ) 



(i) 



where p{f) is the atomic number density field, pg is 
a constant reference density, Sp(f) = p{f) — pi, and 
C 2 (f,f 2 ) = C 2 (\f — f 2 1) j is the two-point direct corre- 
lation function of the fluid, assumed isotropic. 

A general PFC-type functional can be derived from 
Eq. (1), truncated beyond C 2 (\f — f 2 \) : as described in 
Ref. [3J In terms of the rescaled atomic density field 
n(f) = p{f)/pi — 1, an expansion of the logarithm in 
Eq. |l]) generates a functional F = F/(kBTpi) that can 
be written 
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The expansion coefficients w and u are treated as free pa- 
rameters to provide additional model flexibility, and n(f) 
is also typically allowed to assume nonzero average val- 
ues no- Such an expansion is justified for kernels C 2 that 
produce low-amplitude density profiles, which in general 
requires suppression of large wavenumber two-body cor- 
relations. The coefficient w will be set to zero throughout 
this article, as the un 4 (f) /12 term automatically gives 
rise to the effective terms unon 3 (f)/3 and un^n 2 (f) / 
It will also be implied that u = 3 throughout this study. 

Three dynamic equations for n(f) will be considered 
here. The first is a purely diffusive Model B form, 
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where t is dimensionless time. The second equation 
of motion introduces a faster inertial, quasi-phonon dy- 
namic component in addition to diffusive dynamics^, 
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where a and j3 are constants related to sound speed and 
damping rate, respectively. The final equation of motion, 



II. ORDER AND DEFECTS IN PFC MODELS 

In this section, the primary factors controlling defect 
stability in PFC and CDFT models are first outlined 
through an analysis of stacking faults applicable to both 
model types. Three potential solutions to the problem of 
defect instability are examined and shown to be sufficient 
for stabilization of stacking faults in fee crystals. These 
approaches center on multi-peaked correlation functions, 
few-peaked correlation functions with broad effective en- 
velopes, and entropy-driven formulations, respectively. 
Section [TTT] expands on these concepts with a more de- 
tailed examination of specific issues concerning defect 
structure-geometry, interactions, and dynamics. 



A. Large wavenumber or mult i- peaked models 

The efficiency and tractability of PFC models relative 
to CDFT are direct consequences of the central PFC ap- 
proximation; truncation of two-body correlations beyond 
the first few primary correlation peaks in Fourier space 
(which permits the truncated expansion of the logarithm 
in Eq. ([!])). All structural information is consolidated 
into the wavenumber range around the first few primary 
peaks in the structure factor. This still permits control 
of basic structural symmetries, but produces local den- 
sity peaks with broad, sine-wave like profiles rather than 
the sharply peaked Gaussian profiles that emerge from 
CDFT functionals with large wavenumber correlations. 
Though this small wavenumber approximation reduces 
model complexity and increases efficiency by orders of 
magnitude, it can in some cases also reduce defect stabil- 
ity, as will be shown in the present subsection. Methods 
by which defects can be stabilized while still retaining 
the small wavenumber approximation are discussed in 
subsection III Bl 



1. Many-peaked XPFC model 

The effect of large wavenumber correlations on defect 
stability can be demonstrated through an examination 
of the PFC model formulation of Greenwood et alP^Il 
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(XPFC). This approach allows one to stabilize many dif- 
ferent crystal symmetries by constructing a correlation 
kernel with peaks, typically Gaussians, located at the 
first few primary reflections of a given lattice structure. 
The Fourier transformed XPFC kernel can be written 



C 2 [k)i 



' + j ^ e -(^-^) 2 /(2^) e -a 2 fc l 2 /(2p,/3,) ^ 



where i denotes a family of lattice planes at wavenumber 
hi, and a is a temperature parameter. The constants c^, 
Pi, and j3i are the Gaussian width (which sets the elastic 
constants), atomic density, and number of planes, respec- 
tively, associated with the zth family of lattice planes. We 
have introduced the additional constants r and H to per- 
mit C2(k) to assume negative values between peak max- 
ima, which improves numerical stability when used with 
Eq. ([!]). The envelope of all selected Gaussians i shifted 
by the constant — r then composes the final C2(k). For 
example, two such peaks with ratio k<ijk\ = ^4/3 pro- 
duce an fee structure in which the first peak corresponds 
to {111} and the second to {200}, the two primary fee 
reflections. 

The issue of defect stability and its reliance on the 
correlation kernel can be illustrated with an analysis of 
the fee intrinsic stacking fault energy (7isf or SFE) and 
stability as a function of the number of fee reflections 
considered in C 2 (k). The full CDFT of Eq. (I) with 
C 2 (k) given by Eq. ^ was employed, as discussed be- 
low. Faulted crystals were evaluated within a periodic 
simulation cell with (x, y, z) axes in ([112], [110], [111]) di- 
rections, and a single intrinsic stacking fault was created 
by removing one close-packed layer and shifting the lay- 
ers above down by a/y/3, where a = 2\fZi\jk\ — 1.8537 is 
the equilibrium fee lattice constant. The final number of 
close-packed layers Nl must satisfy Nl = 3n — 1, where 
n > 1 is an integer (Nl = 35 was used in most cases 
here). Initial states were relaxed using Eq. In this 
geometry, the stacking fault may be stable, metastable, 
or unstable, with instability unfolding through a shear- 
ing operation r zx or r zy that removes the stacking fault 
and creates a uniformly sheared perfect crystal with total 
shear strain e = (V^A//,) -1 , as shown in Fig. [I] 

Some results are displayed in Fig. [2] along with spe- 
cific parameter values employed. Stacking faults were 
found to be highly unstable in the two-peaked model for 
virtually all meaningful values of (< 2). With five 
peaks, metastable faults were observed, while ten peaks 
produced fully stable faults for which the energy of the 
faulted crystal was very slightly lower than that of the 
sheared state (SF = F$f — ^Shear i$ 0). A 15 peak model 
was also found to be fully stable against shear, with a 
larger negative SF. The relevant structure factors for the 
15 peak model are shown in Fig. [2] These results indicate 
that the stability of stacking faults to shear increases as 
larger-/c correlations are considered. The reason for this 
behavior is discussed later in this section. 

Significantly, the computational efficiency of these 
models decreases rapidly with increasing fc max . The am- 
plitude of the density peaks must be allowed to grow 




(a) 



[111] 

y= [no] 

3= [112] 



(b) 



X 

(c) 



FIG. 1: (Color online) (a) An intrinsic Nl — 35 fee stacking 
fault in a periodic simulation cell; (b) Close-up and full xz 
views of the fault; (c) Full xz view of a sheared/unfaulted 
crystal in the same cell. Yellow lines are guides to the eye. 
These images were generated using Eqs. |2]) and ([9} with pa- 
rameter values no — —0.48, r — —0.63, and B x — 1. 



roughly as /^ ax if the larger-Zc modes are to contribute 
significantly to the free energy of the system. This neces- 
sitates use of the full logarithmic one-body term in the 
free energy, without expansion, since the truncated ex- 
pansion of Eq. p| limits stable amplitudes to relatively 
small values, 0(1). The logarithm greatly reduces the 
maximum allowable time step, while the higher- k modes 
reduce the density peak width (as shown in Fig. [2J top) 
and therefore the maximum allowable grid spacing Ax. 



2. Many-peaked CDFT model 

A similar analysis was attempted using the full RY 
CDFT of Eq. ([l}, with C 2 (k) gi ven b y the Percus-Yevick 
(PY) expression for hard spheres^^ and dynamics given 
by Eq. Our findings indicate that hep crystals in this 
model actually tend to have a slightly lower free energy 
than fee crystals, a fact that apparently has not been 
appreciated before. This leads to a negative 7isf and 
therefore stacking faults that are inherently stable to the 
sheared state. Nonetheless, a perfect edge dislocation in 
a metastable fee crystal was found to split into two well- 
separated Shockley partial dislocations connected by a 
stacking fault, as expected for low 7isf materials (Fig. 
[3|. The reaction is 
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Though both 7isf and 5F in this system are always 
negative or very small and positive, making instability to 
shear virtually impossible, we did find that SF increases 
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FIG. 2: (Color online) Effect of large- k modes on stacking 
fault stability. From top to bottom: Cross-sectional normal- 
ized n(r) maps of the XPFC fee crystal with 2, 5, 10, and 
15 fee reflections from left to right; C2(k) of the PY hard- 
sphere model at pi = 0.9445 and of the XPFC model with 
2, 5, 10, and 15 fee peaks; Log plot of the structure factor 
S(k) = (|n(A;)| 2 ) for faulted and sheared fee crystals in the 15 
peak XPFC model; Linear S(k) plot for faulted and sheared 
fee crystals in the PY hard-sphere model; XPFC 7isf and 
Nl = 35 driving force for instability, 5F — Fsy — i^shear vs. 
fcmax- Symbols in the S(k) plots are Bragg reflection maxima 
of perfect fee crystals in either model. XPFC parameter val- 
ues: pi — 1.0488, cti = 1, a = 0, H = 3, and r = 2. H and 
r were varied from their default values, H = 1 and r = 0, to 
ensure numerical stability of the logarithmic term in Eq. (fTl) . 



as fc max is reduced. This was determined by modifying 
the PY function as C 2 (k) = C% Y (k)e-( k ~ k ™* V 2 for k > 
fc max , while holding Ax constant. The observed trend in 
SF indicates that the underlying link between stability 
and fc max found in the XPFC model is also present here. 



3. Analysis and generalization 

This link can be understood by examining the spec- 
tra of the two relevant states, the faulted crystal and 
the sheared/unfaulted crystal (see Fig. An infinite 
stacking fault can be expressed mathematically as a ID 
step function in the 3D displacement field u{r) of the 
crystal, with the step direction normal to the fault. The 
structure factor of a faulted crystal therefore contains the 




FIG. 3: (Color online) Split fee edge dislocation in the pi = 
0.9445 PY hard-sphere CDFT of Eq. (I). A cutaway view of 
the density field p(r) is shown, with the color scale truncated 
above p(f) = 1, though p(f) ~ 7800 at the peak maxima. 
Left: Diagrammatic representation of the dissociation of a 
perfect dislocation b into Shockley partials bi and b^. 

same 1/k 2 line-shapes characteristic of scattering from a 
surface. In this case each lattice reflection fcW, except 
the set of {111} planes parallel to the stacking fault, will 

possess its own set of l/(k z — k^) 2 modes, and these 
modes will only extend in the reciprocal space direction 
perpendicular to the plane of the fault, k z . The spher- 
ically averaged structure factor S(k) thus exhibits a se- 
quence of broadened fee reflections, each with effective 
l/(k — k^) 2 line-shapes of roughly the same width 30 . 

The structure factor of a sheared perfect crystal also 
exhibits a form of peak broadening, but in this case the 
degree of broadening is proportional to k. Depending on 
the orientation of shear strain e, or deformation in gen- 
eral, the spacing d between certain lattice planes within 
any given family changes to some value d! . A simple 
analysis shows that the resulting shift in /c-space under 
shear is 

( 1 1 \ 1 — COSCJ 7 , v 

5k = 2ir - - - = k 8 

\ d' d J cos uj 

where uj = arctane. Thus the primary \ow-k reflections 
remain sharp and/or slightly split, while the degree of 
shift or effective broadening increases in proportion to k 
due to this growing decoherence effect across high-index, 
\ow-d planes. 

The driving force for stacking fault instability is pri- 
marily a function of SF = Fsf — ^Shear- When a low- 
mode kernel is used, Fshear will tend to be smaller than 
Fsf because the slightly shifted low-k modes of a sheared 
crystal induce a small energy cost relative to that of the 
long-range line-shape tails of a faulted crystal, provided 
that the shear strain is not too large. As k max increases, 
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this situation eventually reverses. The high order reflec- 
tions from the sheared lattice become broader or less 
coherent than those of the faulted crystal such that if 
enough high-/c modes are considered, stacking fault sta- 
bility can eventually be attained when Fshear > ^sf- 
This effect can be seen in the data of Fig. [2] We em- 
phasize that the high-/c modes do not necessarily decrease 
7isf or Fsf, but that the increase in Fshear alone leads to 
a dramatic increase stacking fault stability. In the geome- 
try of our simulations, the critical wavenumber at which 
^sf = ^Shear increases with Nl since e = {\/2Nl)~ 1 1 
but for fixed Nl stability increases with fc max . The exact 
wavenumber at which stability is attained (/c max ~ 18 
in this case) also depends on model parameters, but the 
basic trends should generally hold. 

Similar arguments can be constructed for other types 
of defects by reconsidering for each case the two generic 
parameters controlling stability, driving force and energy 
barrier. The driving force for instability is controlled by 
both the energy of the defected state and the energies of 
the states resulting from defect removal. In the case just 
examined, this force was quantified by 5F, the energy 
difference between the faulted state and the competing 
sheared state. The magnitude of this force is therefore 
deeply tied to the precise form of (^(fe) and the energy 
cost that it imposes on defect line-shapes, a few examples 
of which are shown in Fig. [4] These spectra were obtained 
using the modulated version of the original PFC model 
presented in Section III A[ though the basic features do 
not depend on the model used. The four crystals shown 
have roughly equivalent defect densities and their struc- 
ture factors exhibit similar line-shapes, with the stacking 
fault producing the greatest degree of broadening. It is 
apparent that in general, broad correlation kernels will 
reduce the energies of all of these defect structures and 
therefore the driving force for their instability, as demon- 
strated explicitly in the next subsection. Alternatively, 
stacking faults were stabilized in the previous example by 
instead increasing the energy of the post-defect-removal 
sheared state through high-/c correlation effects. This 
simply illustrates that the driving force can be controlled 
by varying the energy of the initial state, the final state, 
or both through the form of C^W- 

The second parameter relevant to stability, the energy 
barrier for defect removal or topological protection, is 
controlled by the nature of the process by which a given 
defect can be removed. The shearing operation that re- 
moves a stacking fault requires only small local trans- 
lations of the density peaks. No peaks are created or 
removed, thus no long-range mass transport is involved, 
and the resulting energy barrier is small. Stability in 
such a case therefore requires that the driving force be 
minimized, since the energy barrier is in general more dif- 
ficult to manipulate. Removal of dislocations and grain 
boundaries on the other hand requires more elaborate 
transformations that often involve climb, with redistri- 
bution of atoms by long-range diffusive mass transport 
(adding or removing density peaks in the case of a PFC 



Perfect fee, H =0.00 
Intrinsic Stacking Fault, H =0.00 - 
Edge Dislocation Pair, H =0.08 - 
Screw Dislocation Quadrupole, H =0.08 - 
Grain Boundary Pair, H =0.00 




FIG. 4: (Color online) Structure factors of fee crystals with 
various types of defects. Note that the most compact struc- 
ture, the planar stacking fault, produces the overall broadest 
line-shapes. These spectra were generated using Eqs. (pi and 
(12} with parameter values n = -0.48, ~ - n ao ° r - 1 
ao = 1/2, and ko = 6.2653. 



-0.63, B x 



1, 



model). These defects thus have greater inherent sta- 
bility than stacking faults due to both their somewhat 
narrower line-shapes and their generally larger energy 
barriers for removal. 

The results presented in this section demonstrate that 
certain defect structures can be stabilized in PFC models 
by building structural information back into the model 
with additional correlation peaks. The resulting multi- 
peaked models require the full logarithmic free energy 
functional and therefore suffer from marked inefficiencies 
and greater complexity relative to existing PFC models 
studied in the literature. However, the sources of defect 
instability uncovered here will prove useful in analyzing 
and constructing few-peaked PFC models that both re- 
tain efficiency and stabilize relevant defects. 



B. Small wavenumber or low-mode models 

The inefficiencies associated with large wavenum- 
ber models can be avoided by maintaining the small 
wavenumber PFC approximation and giving closer treat- 
ment to the details of the one or few correlation peaks 
employed. We have successfully stabilized stacking faults 
using four variations of PFC with the small wavenumber 
approximation, though in general this requires some de- 
gree of parameter tuning and/or slight modification of 
C 2 (k). 

Perhaps the simplest way to improve defect stability in 
this approximation is to use a broad C^CO in the region 
of the first few primary lattice reflections. This in gen- 
eral lowers the energy associated with defect line-shapes 
while having little to no effect on the energy of, for ex- 
ample, uniformly sheared states or undeformed crystals. 
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A sufficiently broad envelope can thus stabilize the 1/k 2 
stacking fault modes relative to the first few sheared state 
Sk shifts. The main drawback is that the stability of the 
equilibrium crystal structure tends to decrease as C 2 (k) is 
broadened. For example, a two-peaked XPFC fee model 
can prefer hep, bec, rod, or lamellar structures over fee 
for sufficiently large c^. Thus defect stability is not given 
freely, it comes inversely bound to crystal stability. This 
trade-off must be managed by ensuring that the primary 
reflections of the desired crystal structure continue to be 
preferred over those of competing symmetries that may 
exist within a broad C 2 (k) envelope. Specific parame- 
ter values related to the width of C 2 (k) (elastic moduli) 
therefore become restricted, but a workable balance be- 
tween defect stability and crystal stability is still often 
attainable in the small wave number PFC approximation, 
as demonstrated in the following subsections. 



1. 



PFC model 



We begin with the original PFC model of Elder et 
alP^, for which the correlation kernel and the primary 
fee and hep reflections are shown in Fig. [5] This PFC or 
Brazovskii 32 functional specifies a kernel 



2\2 



C 2 {k) = -r+l-B x (l-k 2 ) 



(9) 



where B x is a constant proportional to the solid-phase 
elastic moduli and k = fc/(27r). This C 2 (k) produces 
equil ibrium fee structures within a certain parameter 
ran g^33lMl 5 b u t the energy of the hep crystal is neces- 
sarily always very close to that of fee, and is sometimes 
lower. This means that the fee 7isf is always vanish- 
ingly small in the original PFC model, which naturally 
stabilizes stacking faults relative to shear. Perfect dislo- 
cations, as with the CDFT hard-sphere model, split into 
two well- separated Shockley partials and a stacking fault, 
as shown in Fig. [6j Note that, for display purposes, the 
defects shown in Fig. [6] were generated using a modified 
version of Eq. ([9|, discussed in Section [ill A The small 
value of 7isf produced by Eq. (|9| leads to an equilib- 
rium splitting distance ~ 32a for b — a/2 [110] edge dis- 
locations, which is significantly larger than values < 10a 
observed in typical single component fee metals. 

The broad C 2 {k) function of the original PFC model 
therefore naturally stabilizes defects, but is relatively lim- 
ited in terms of flexibility in controlling crystal symme- 
try and elastic properties. One is restricted to the deep 
quench region of fee stability, where small changes to the 
elasticity parameter B x tend to destabilize fee, and the 
inherent 7isf is extremely small. It is possible to tune 
7isf within a small range by varying the envelope width 
B x , but a more effective approach is to explicitly enhance 
or suppress specific hep and fee modes, as demonstrated 
in Section HTT1 



Original PFC — 




FIG. 5: (Color online) Sample PFC kernels that produce 
stable stacking faults, with relevant fee and hep reflections 
also indicated. Inset: Close-up view of the peak maxima. 
Original PFC parameters: n = -0.48, r = -0.63 and B x = 
1. XPFC* parameters, Eq. (To): n = -3/10, r = -9/40, 
W2 = 1/50, w A = 49/50, an ~ 1, a = 0, and H ~ 1.0625. Wu 
et al. parameters: no = —3/10, r = —9/40, Ri = 1/20, and 
B x = 1. Note that for display purposes all kernels have been 
shifted by —Suq to maintain consistency with Eq. |l]), and 
that k has been rescaled in the XPFC* and Wu et al. kernels 
to match the equilibrium fee reflections of the original PFC 
kernel. 



2. Low-mode XPFC model 

XPFC models with only a few low-A: peaks can also be 
tuned to support less protected defect structures, though 
some modifications seem to be necessary in the case of 
fee stacking faults. As is increased above ~ 1 or 2, 
Gaussian kernel peaks generally become too broad to 
maintain a consistent crystal structure before they are 
broad enough near their maxima to sufficiently support 
the 1/k 2 stacking fault modes. We have obtained better 
results in this regard using a modified Gaussian shape 
function, 



C 2 (k) z = 



He 



(10) 



where the coefficients w 2 + W4 = 1 set the relative weight 
of k 2 versus k A shape functions. This formulation still 
permits control of the phase diagram in the original spirit 
of the model as well as the magnitude and anisotropy of 
elastic constants, though the range of parameter values 
sufficient for stability remains restricted. Values in the 
vicinity of w 2 = 1/50, W4 = 49/50, = 1, a = 0, 
H — 1, r = —9/40, and no = —3/10 have been found 
sufficient to stabilize stacking faults in a three-peaked 
fee model with two fee reflections at k 2 = ^4/3fci and 
one commensurate hep reflection at k% = ^/41/36fci (see 
Fig. [5|. Similar results have been obtained using only 
the two fee reflections with their maxima connected by 
any nearly linear bridge function. 
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FIG. 6: (Color online) Dissociated b = a/2 [110] edge and 
screw dislocat ion s in the (modulated) original PFC model of 
Eqs. p]) and ([T^. The lower images show cutaway views of 
n(r), illustrating the mixed Shockley partials connected by a 
single intrinsic stacking fault. The upper images show atomic 
representations of the peak positions in n(f): red atoms have 
local hep coordination (stacking fault), gray atoms have local 
noncrystalline coordination (defect cores), and atoms with fee 
coordination are not shown. Parameter values no = —0.48, 
r -0.63, B x = 1, a = 1/2, and k 6.2653 were used. 



where Qi = \/4/3 is the wavenumber of the second mode 
and Ri sets its height. We find that quenches of interme- 
diate depth tend to produce fee crystals with moderate 
7isf 5 and that 7isf decreases with increasing |no|. This 
is because the energy of hep approaches and may become 
lower than that of fee as the quench becomes deeper. As 
an example, the parameter set r = —9/40, R\ = 1/20, 
B x = 1, and no = —1/4 produces metastable stacking 
faults for Nl = 35 and perfect dislocations that do not 
split into partials. As no is lowered, Nl = 35 stacking 
faults eventually become fully stable, with no = —3/10 
producing a very low 7isf and widely split partial dis- 
locations (see Fig. [5|. 7isf also varies rapidly with B x , 
again permitting only a small degree of elastic moduli 
tuning. 

The source of defect stability in this PFC model is 
apparent from examination of the kernel around the pri- 
mary fee modes. Elasticity parameters B x < 1 produce 
broad profiles with very shallow wells between to the two 
modes, while B x > 1 produces narrowing profiles with 
increasingly deeper wells. Using the parameter values 
above at no = —3/10 and B x = 2, the stacking fault line- 
shapes are already unfavorable enough to fully destabilize 
Nl = 35 faulted crystals. This two-mode model comes 
with similar restrictions in terms of acceptable param- 
eter ranges to those of the original and XPFC models, 
but does provide greater control of fee stability than the 
original PFC model, allowing use of shallower quenches 
closer to the linear elastic regime. 



These findings highlight the conflict between crystal 
stability and defect stability in small wavenumber mod- 
els. Very narrow crystal-stabilizing XPFC kernels desta- 
bilize defects of all but the most protected types, some- 
what broader kernels still tend to destabilize suscepti- 
ble planar fault structures, while very broad fully defect- 
stabilizing XPFC kernels ultimately resemble the simpler 
kernels of the original PFC model. These inherent trade- 
offs suggest a limitation to the range of atomic-level crys- 
tal and defect structures that low-mode, two-body PFC 
models can capture, though the ultimate limit appar- 
ently lies at some level of complexity beyond that of the 
primary defect structures in fee crystals. 



S. Wu et al fee PFC model 

The two-mode fee model introduced by Wu et al.^also 
produces stable stacking faults within certain parameter 
ranges. To facilitate tuning of 7isf and stability, we find 
it helpful to retain a bulk modulus coefficient B x analo- 
gous to that of the original PFC kernel, 



c 2 (*o 



2\2 



l-B x (l-k 2 ) 



B x 



(11) 



4- Vacancy PFC model 

The vacancy PFC or VPFC model typically employs 
the same C^fc) as the original PFC model, but adds a 
strong local nonlinear cutoff term that destroys the in- 
herent one-mode nature of inhomogeneous density states 
and produces relatively autonomous indi vidu al density 
peaks separated by regions of zero 

densit} |3613Z| Thig 

can 

lead to stabilization of individual lattice vacancies and 
in general produces a Brownian hard-sphere-like system 
in which fee and hep structures have lowest free energy, 
with a slight preference toward fee. When vacancies are 
stabilized in this way, packing effects become more im- 
portant and the natural similarities between fee and hep 
in this regard should lead to low 7isf systems. We find 
that a low 7isf, comparable to that of the original PFC 
kernel, is indeed obtained with the VPFC description, 
though in this case we attribute the low value primarily 
to entropic or packing effects. In this sense, the VPFC 
model is likely the most defect-friendly of existing PFC 
models, but one must reinterpret the time scales accessi- 
ble to this description since the explicit vacancy diffusion 
mechanism has been partially reintroduced. Nonetheless, 
stacking faults are quite stable in this formulation and 
perfect dislocations split properly into partials joined by 
a stacking fault. 



8 



III. SURVEY OF FCC DEFECT PROPERTIES 
IN A SMALL WAVENUMBER PFC MODEL 

The findings presented in Section [TT] suggest that many 
types of crystalline defects observed in real materials, 
even those with relatively low inherent stability, can be 
stabilized within existing PFC models. The general re- 
quirement for stability under the small wavenumber ap- 
proximation involves a sufficiently broad correlation ker- 
nel maximum, which leads to restrictions on the parame- 
ter ranges that provide both defect and crystal stability. 
Nonetheless, these restrictions appear to be manageable 
in the case of close-packed crystals, and closer examina- 
tions of fee systems reveal that a wide variety of realistic 
defect properties naturally emerge within these parame- 
ter ranges. Some of these properties are discussed in the 
present section. 

Emphasis is placed on fee crystals and tuning of the 
stacking fault energy to generate either undissociated 
perfect dislocations (high 7isf) or dissociated Shockley 
partial dislocations bound by stacking faults (low 7isf)- 
Various fundamental properties of both edge and screw 
dislocations in these configurations are examined. 



A. Modulation of the original PFC model 




FIG. 7: (Color online) Modulated original PFC kernel of Eq. 
(T2]) where n -0.48, r -0.63, B x = 1, a = 1/2, and 
ko = 6.2653. Ho = 0, 0.04, and 0.08 are shown. Inset: 7isf 
obtained with the given modulated kernel vs. Fh cp — F{ cc and 
Ho. 71 corresponds to Nl = 53 and 72 to Nl =35. As in 
Fig. [5] all kernels have been shifted by — 3no- 



stable or metastable for Ho < 0.05, and very little system 
size dependence is observed in 7isf for Nl > 35. 



As discussed in Section II B 1 , the dissociation width of 



a perfect dislocation in the original PFC model is larger 
than that of most single component fee metals, due to 
the small value of 7isf relative to the fee elastic moduli. 
This behavior can be tuned by various means. Since 7isf 
is strongly linked to the difference in free energy between 
commensurate fee and hep states, one method to control 
7isf is to selectively promote or suppress the primary 
hep reflections that do not coincide with those of the fee 
crystal. For example, the commensurate hep reflection 
at k 2tt ^41/12/ a can be suppressed by subtracting a 
small XPFC-type Gaussian centered at this wavenumber 
from the of Eq. ([9]), such that the two fee reflections 

are not affected, 

C 2 {k) = -r + l-^(l-P) 2 -i^ e" (/c " /co)2/(2a ° ) . (12) 
The constants Ho and ao are analogous to those of Eq. 

The original PFC kernel modified in this way is shown 
in Fig. [7| The result is greater stabilization of fee due 
to the increase in Fh cp , and a corresponding increase in 
7isf as the height of the subtracted Gaussian Ho be- 
comes larger. This C^fc) is similar to that of Wu et al. 
with variable B x , though some subtle differences are ap- 
parent. Only 7isf of the equilibrium fee state changes 
as Ho is varied, producing an essentially fixed reference 
system with tunable 7isf- The standard phase diagram 
is of course modified for nonzero Ho, proportionally ex- 
panding the region of fee stability. The dependence of 
7isf on Hq is also shown in Fig. [7J Nl = 35 faults are 



B. Dissociation width vs. 71SF 

This method of controlling 7isf can be used to exam- 
ine defect properties that depend on stacking fault en- 
ergy, and to compare functional dependencies with those 
predicted by continuum elastic theories. For example, 
the elastic strain energy of a fee crystal containing a per- 
fect dislocation with b = a/y/2 is lowered by dissociation 
into two Shockley partials with b p = aj\f§ and a stack- 
ing fault, as described by Eq. (|7|). This is confirmed by 
the Frank criterion, b 2 = a 2 /2 > 2b 2 = a 2 /3. Contin- 
uum elastic theories provide predictions for the equilib- 
rium separation d eq between the two resultant Shockley 
partials. The long-range elastic energies of the partials 
produce a repulsive force per length of dislocation line 
proportional to while the energy cost of the fault in- 
creases with d, producing an attractive force per length 
proportional to 7isf- The separation that balances these 
two forces for an elastically isotropic material is given 
b^ 



eq 



1 



2v cos 2/3 \ /xfep 



87T7ISF 



(13) 



where v is the Poisson ratio, j3 is the angle between the 
Burgers vector of the perfect dislocation and its line di- 
rection, and /i is the isotropic shear modulus of the crys- 
tal. Similar though more complicated results can be de- 
rived using anisotropic elasticity or various approxima- 
tions thereof^. 
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FIG. 8: (Color online) Equilibrium splitting distance 
between Shockley partials within dissociated perfect b = 
a/2 [110] edge and screw dislocations, shown as a function 
of 7j^p. Lines represent the predictions of isotropic and 
anisotropic continuum elastic theories, with no adjustable pa- 
rameters. Model parameters are the same as those of Fig. [7] 



We have measured d eq for both edge and screw dislo- 
cations in the modulated PFC model, as a function of 
7isf- Base parameter values no = —0.48, r = —0.63, 
and B x = 1 were used, giving an fee lattice constant 
a = 1.8537. The crystal orientation of Section [TT] was 
employed (see Fig. [6|, though in this case with peri- 
odicity broken in the z-direction by a thin layer of liq- 
uid. A single perfect b = a/2 [110] edge dislocation with 
line direction [112] was initialized at (L y /2,L z /2) using 
system dimensions (L Xl L yi L z ) = (16,3986,1820) and 
Ax = ^3/512 (203,668 atoms). Perfect b = a/2[110] 
screw dislocation dipoles with (110) line directions were 
initialized similarly at (L x /^L z /2) and (3L x /4, L z /2) 
with (L x ,L y ,L z ) = (4278,10,1820) and Ax = a/V^00 
(107,692 atoms). The symmetry of the screw dislocation 
displacement field in the ^-direction necessitates a dipole 
configuration. In both cases the diffusive dynamics of 
Eq. (§ were used. Results are shown in Fig. [8j including 
comparisons with Eq. (13) and the available anisotropic 
theories. 

The expected linear trend is clearly observed for both 
dislocation types, and the agreement with linear elas- 
tic predictions is excellent. The fully isotropic the- 
ory very slightly overestimates d eq , while the isotropic 



'eq? 

theories with approximate anisotropic moduli of Voigt 
and Reuss^J and the fully anisotropic theory (which 
is tractable only for the screw dislocation) match the 
data extremely welPn We note that it was necessary 
to simulate relatively large systems to eliminate signif- 
icant boundary and image effects. The periodic image 
forces nonetheless induce a small bias toward larger d 
that grows with d and appears to become appreciable 
for d > 25a at this system size. Otherwise, the largest 
proportional disagreement is seen at small values of d 



where the linear elastic predictions are likely somewhat 
inaccurate due to core overlap effects. We believe that 
these results provide an important validation of funda- 
mental defect energetics and interactions within the PFC 
description. 



C. Dissociation width vs. applied (Escaig) strain 

When the glide-inducing shear strain e zy is applied to 
either of the dissociated dislocation lines shown in Fig. 
|6j the two partials will glide together with a nearly con- 
stant separation (except at large velocities as discussed 
below). This is because the components of their respec- 
tive Burgers vectors in the shear direction [110] have the 
same sense. This is not the case for shear strain e zx . The 
two partials have components in this shear direction [112] 
with opposite sense, which produces glide forces in oppo- 
site directions that extend the faulted region in between. 
This particular strain orientation is sometimes referred 
to as the Escaig strain. The new equilibrium separa- 
tion for fixed stress r zx can again be calculated within 
isotropic elasticity theory, and for a general dislocation 
one obtaimP 



d 



tt(27] 



ISF 



6„ | sin 6 2 



sin Vi | 



(14) 



where f (61,62) — cos 61 cos 62 + sin 61 sin 62 / ( 1 — v) , and 
#1, 62 are the angles between the overall dislocation line 
direction and the Burgers vector directions of the leading 
and trailing partials, respectively. For the case of a screw 
dislocation, where 1 = -30° and <9 2 = 30°, Eq. ^ 
reduces to 



d = 



2-Zv 



\ib 2 v 



8tt(1 - v) 7isf - r zx b p /2' 



(15) 



Shear strain e zx was applied to the screw dipole config- 
urations described in the previous subsection by adding a 
penalty function to the first few crystalline surface layers 
near the liquid boundaries and translating the penalty 
field at some constant rate e zx to drive the external de- 
formation. An affine shear deformation was also applied 
to the entire crystal at each time step, and the inertial 
equation of motion, Eq. Q, with a = 1 and (3 = 1/100 
was used. These features together should ensure that the 
resolved shear strain at the dislocation remains as close 
to the applied shear strain as possible. Selected fixed val- 
ues of e zx were held periodically to allow the system to 
fully relax to a steady configuration. Though a reliable 
method of directly quantifying stress in our simulations 
is currently lacking, the assumption of linear elasticity 
will produce accurate stress-strain conversions whenever 
plastic flow / stress dissipation is negligible and strain is 
not exceedingly large. We can confidently apply this as- 
sumption here (r zx = C zx e zx ) to obtain results in terms 
of stress, which can be directly compared with Eq. (15). 
Expected errors are < 1 — 2% for e zx < 0.05. 
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FIG. 9: Divergence of the equilibrium splitting distance 
between Shockley partials within dissociated perfect b = 
a/ 2 [110] screw dislocations, shown as a function of e zx , the Es- 
caig shear strain. Lines represe nt t he predictions of isotropic 
continuum elastic theory, Eq. (15), as discussed in the text. 
Model parameters are the same as those of Fig. [7] 



Simulation results are compared with the predictions 
of Eq. (15) in Fig. [9] The lines show Eq. (15) with elastic 
parameters as determined from simulations. The agree- 
ment is in general good for all values of i^o, with the ex- 
pected trend of diverging separation clearly visible. The 
strain at which divergence appears to occur agrees well 
with the predictions for small 7isf> though not surpris- 
ingly shows increasing deviation as 7isf and e zy become 
large. These findings demonstrate that non-trivial defect 
properties associated with interactions between defects 
and applied stresses can be accurately modeled with this 
approach. 



D. Peierls stress/strain for dislocation glide 

Dislocation motion is the dominant microscopic ele- 
ment of plastic deformation in most crystalline and poly- 
crystalline metals. The minimum strain required to move 
a dislocation, the Peierls strain ep, is therefore a funda- 
mental material property that can qualitatively alter a 
material's macroscopic plastic response, yield stress, etc. 
fee metals, in which dissociated dislocations and stacking 
faults are prevalent, typically exhibit a characteristic slip 
activity that derives directly from the nature of the dis- 
sociated dislocation structure and the anisotropy of its 
Peierls barrier. Stacking faults are stable only on {111} 
planes, and the Peierls barrier for motion of the partials is 
normally very small or negligible within the same {111} 
slip planes. The result is that a large majority of slip 
activity occurs only within {111} planes, and a regime 
of qualitatively different plastic response emerges when 
cross-slip between {111} planes becomes active. Thus 
a faithful description of the underlying Peierls barrier is 



central to any model of crystal plasticity, and will nat- 
urally give rise to many of the secondary features and 
processes that emerge from the fundamental mechanism. 

Accurate measures of Peierls stresses rp in fee metals 
have proven difficult to obtain through atomistic com- 
putations, due to the smallness of rp relative to typical 
boundary image stresses in finite size simulations 41 . Typ- 
ical estimates thus vary greatly in the literature, but the 
most reliable values measured from molecular dynam- 
ics (MD) simulations are believed to be on the order of 
10 _5 -10 _4 /i for fee e dge disl ocations and 10 _4 -10 _3 /i for 
fee screw dislocation d 38 * 41 * 4 ^. Experimentally determined 
values also vary greatly, from ~ 10 -6 — 10 _3 /i, depending 
on the method employed 4 ^ 4 ^. Nonetheless, within these 
general ranges, screw lines are expected to have larger 
barriers than edge lines, and unsplit perfect dislocations 
to have larger barriers than split partial configurations. 

The Peierls strains for glide of edge and screw disloca- 
tions, dissociated and undissociated, were measured by 
applying shear to each system as described in the previ- 
ous subsection. Eq. Q with a = 1 and j3 = 1/100 was 
again employed, though in this case shear strain e zy was 
applied. This orientation, unlike e ZXl produces uniform 
glide of both dislocation types, with split partials gliding 
in the same direction, ep was defined as the magnitude 
of applied strain at the instant a given defect has glided 
a distance equal to its Burgers vector magnitude b. 

The simulation results shown in Fig. [To] bear out all of 
the expectations noted above with only relatively minor 
quantitative deviations. All dislocations exhibit the same 
roughly \fl barrier strain dependence for sufficiently large 
shear rates. The cause of this behavior is discussed be- 
low, but our primary interest is in the limiting values 
obtained at small shear rates. Since the amount of dis- 
location motion in these simulations is negligible, we 
can again apply linear stress-strain relations to report 
Peierls stresses rather than strains (rp ~ /iep). The 
unsplit screw Peierls stress approaches a value within 
the expected range rp ~ 6 x 10 _4 /i, while the aver- 
age Peierls stress for the split screw is slightly lower, 
rp ~ 3.5 x 10 _4 /i. The measured values for the leading 
and trailing screw partials differ by a small amount. This 
may be a consequence of weakly asymmetric periodic im- 
age forces generated by the simulation boundaries, or of 
an initial d value that is incommensurate with the peri- 
odicity of the Peierls potential, as discussed in Ref. [45] 

The widely split edge dislocation exhibits much less 
deviation between leading and trailing rp values, though 
the apparent limiting value is not smaller than that of 
the perfect edge dislocation to within expected error. 
Both configurations approach a value that appears to 
extrapolate to rp ~ 3 x 10 _5 /i, also within the expected 
range for typical fee edge dislocations. It was observed 
that edge dislocations tend to glide with a relatively uni- 
form, continuous progression, while screws exhibit a more 
pronounced stick-slip hopping motion with each unit of 
translation, indicative of a larger Peierls barrier. These 
results together demonstrate that the fundamental fea- 
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FIG. 10: Measured Peierls strain for glide of dissociated 
and undissociated edge and screw dislocations, shown as a 
function of applied shear strain rate. The solid and dashed 
lines, which are indistinguishable, are fits to the functional 
forms presented in the text. Model parameters are the same 
as those of Fig. [7] 



tures of fee Peierls stresses are captured quite reasonably 
by this relatively simple PFC model. 

The apparent relation ep ~ \f e for larger shear rates 
can be obtained as follows. Based on our results, there 
appears to be an intrinsic time scale r associated with 
relaxation of an entire defect structure, inluding the core 
region and the long-range displacement fields. If shear 
is applied slowly relative to r, the core region and the 
long-range displacements can evolve cooperatively such 
that the inherent Peierls barrier of the given dislocated 
system, 5Fp, is realized. If e _1 is large relative to r, 
then the short-range and long-range components can- 
not relax with optimal synchronization. This introduces 
nonequilibrium effects and alternate relaxation processes 
that lead to measurement of some larger effective barrier 
5 Fes = SFp + <^neq- 

In the PFC formulation, slow relaxations near equilib- 
rium are diffusive whether one uses Eq. (|3) or Eq. @. 
Thus the intrinsic defect time scale r can be approxi- 
mated as an exponential relaxation time, and neglecting 
other relaxations, Fneq should be, to a first approxima- 
tion, an exponential function of the dimensionless pa- 
rameter re. We therefore may write SF^eq ~ e re — 1 or 
SFneq ~ r £ m the small re limit. Combining this rela- 
tion with a linear elastic response, SF E ^ ~ e 2 , we obtain 
ep ~ \J SpEf? ~ V SFp + ere where c is a constant. This 
form and the corresponding unexpanded form both fit 
our results well, indicating that the competition between 
rate of applied shear and rate of diffusive strain field re- 
laxation leads to an increased effective Peierls strain with 
relevant control parameter re. The observed system size 
dependence of ep is also consistent with this expression. 



FIG. 11: (Color online) Measured dissociation width of a 
gliding edge dislocation vs. dislocation velocity for various 
applied shear strain rates. Steady-state widths and velocities 
at each shear rate are separately indicated. The solid line is 
a fit to the functional form of the drafting theory presented 
in the text, and the dashed line is an illustrative example 
of the behavior expected from semi-relativistic effects (d ~ 
y^l — v 2 /a 2 ). Model parameters are the same as those of 

Fig. [3 



E. Dynamic contraction of dissociated dislocations 

When the glide-inducing stress r zy is applied to a dis- 
sociated dislocation, the Peierls barrier is overcome quite 
quickly and glide is initiated within the slip system of 
the defect. The structure of the extended defect may 
then be altered by dynamic effects associated with navi- 
gation of the Peierls potential as well as drag or damping 
forces of various origins. Employing the dynamics of Eq. 
Q, we applied fixed shear strain rates to the dissoci- 
ated edge and screw dislocation configurations already 
described, and allowed each to approach a steady-state 
glide velocity v, while montoring the separation between 
paired partials d. Similar results were obtained for both 
dislocation types, and those for the Ho = edge dislo- 
cation are shown in Fig. [TT] For relatively low steady- 
state velocities (v < 0.1a), d was found to exhibit regu- 
lar oscillations in time. These oscillations, which can be 
seen in the lowest shear rate data displayed in Fig. [TTJ 
are not unlike t he so- called breathing modes observed in 
MD simulation d 46 * 47 *, though their period appears to be 
significantly longer and more regular in these PFC sim- 
ulations. 

At higher velocities a dramatic decrease in d was ob- 
served, followed by a return approximately to the initial 
d upon release of strain. Since the shear strain e zy does 
not significantly alter 7isf or the repulsive force between 
static paired partials, this contraction appears to be a 
fully dynamic effect. We believe that the primary mech- 
anism at work is a difference in the frictional drag force 
experienced by each partial at high PFC velocities or, 
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similarly, a reduction of the drag force exerted on the 
extended defect structure as a whole. Either scenario 
may result from the difference in screw-edge character 
between the two partials 38 or may be driven by a co- 
operative drag-reduction process analogous to the phe- 
nomenon of drafting or slipstreaming in fluid dynamics. 

At low velocities, where the effect of drag is small, v 
is determined primarily by the height of successive bar- 
riers in the Peierls potential. Thus any forces that might 
change the dissociation width are small, though they may 
become large enough to generate low-amplitude breath- 
ing modes. As v increases, the Peierls potential becomes 
less relevant and the effect of frictional drag becomes 
large. In this regime, a dissociated defect may experi- 
ence a larger total drag force than an undissociated one 
due to its greater extension. The total drag force may be 
reduced, thus increasing v at a given strain, by reducing d 
such that the trailing partial assumes an optimized lower- 
drag position within the leading wake or slipstream. At 
large enough velocities, the reduction in drag force expe- 
rienced by the full defect under contraction may compen- 
sate for the increase in local defect energy. The energy 
of the entire system is reduced by the excess strain relief 
that is a consequence of faster dislocation glide. 

We can modify the analysis given in Ref. 38 to incorpo- 
rate such an effect into the predicted steady-state value of 
d. The steady-state force balance between gliding paired 
partials can be written^ 



7isf J rB = D 1 
7isf + D 2 = C 



A 

~d 
A 

~d 



(16) 



where A is an abbreviation of the terms multiplying 7 IS p 
in Eq. (13), B and C are the forces per unit length on b\ 



and 62, respectively, due to the applied shear, and D\ and 
D2 are the damping forces per unit length on the moving 
partials. First we assume that both partials experience 
the same force due to the applied shear strain (B — C) 
and that the nominal drag force on each PFC disloca- 
tion obeys a Stokes' law form [Di ~ Cdv), as shown in 
Ref. [9j The drag coefficient should roughly be given by 
Cd — fi/M, where M = v/e ~ 3 is the partial dislocation 
mobility. Finally, since we assume that D2 decreases with 
d (D2 = Cdvf(d)) a plausible functional form for f(d) 
must be proposed. We will use f(d) = Sd/ (Sd + c v v) 
where Sd = d + do, and do, c v are constants. 

A fit to the resulting equation for d vs. v is shown in 
Fig. [TT] The only adjustable parameters in the fit are do 
and c Vl and the agreement is very good. Semi-relativistic 
contraction of strain fields in the glide direction has also 
been predicted for dislocations moving at velocities near 
the sound speed a of the crystal. A curve illustrating 



this effect is shown in Fig. 11 , where it is assumed that 
both partials experience the same contraction effect and 
reduce d accordingly. This mechanism does not have a 
large effect until the sound speed is approached and its 
form is not consistent with our results. Furthermore, we 



observe similar behavior to that shown in Fig. 11 when 
dynamics are given by Eq. (|3|, which does not introduce 
a sound speed. 

The PFC description is generally best suited for exam- 
ining behavior under relatively low strain rates or driv- 
ing forces, where dislocation velocities tend not to closely 
approach the sound speed. Nonetheless, a better under- 
standing of the nature of the interaction between moving 
dislocations and the quasi-phonons described by Eq. Q 
would be useful in terms of confirming that the artifi- 
cially low sound speeds employed in PFC simulations do 
not qualitatively alter the low and intermediate velocity 
dislocation dynamics. Details of these issues are deferred 
to a future publication. 



IV. DISCUSSION AND CONCLUSIONS 

The primary factors controlling defect stability in PFC 
models have been examined, and it has been demon- 
strated that broad correlation kernels or elastically soft 
crystals produce the greatest defect stability. Maximally 
broadened kernels appear to be necessary for stabiliza- 
tion of certain defects such as planar faults in the small 
wavenumber PFC approximation. Higher wavenumber 
correlations with narrower peaks were also shown to im- 
prove stability in some cases, but this feature leads to 
greatly reduced model efficiency. The inherent conflict 
between crystal stability and defect stability in broad 
kernel models suggests potential limitations to the com- 
plexity of structures that can be described by such mod- 
els. Any such limitations were shown to be non-factors in 
the case of the primary defect structures in close-packed 
fee crystals. 

Stacking faults with low inherent stability were stabi- 
lized in four PFC variants, indicating that considerable 
defect stability can be obtained without losing crystal 
stability. The central defect structure in fee plasticity, 
dissociated partial dislocations with stacking faults, has 
been examined in some detail and shown to be well- 
described by PFC methods. The dependence of the 
equilibrium dissociation width on stacking fault energy 
has been shown to agree with continuum elastic predic- 
tions for both edge and screw dislocations, under zero 
and nonzero applied external stresses. Peierls stresses 
for glide of both dislocation types have also been mea- 
sured and shown to fall within the typical ranges and rel- 
ative magnitudes determined from experiments and MD 
simulations. Contraction of gliding pairs of partials has 
been observed at high velocities and argued to be a con- 
sequence of large frictional drag forces on widely split 
partials. 

The findings presented in this article are intended to 
lay the groundwork for larger-scale PFC studies of plas- 
ticity in fee crystals, with the potential to examine exper- 
imentally relevant strain rates, nonconservative/climb- 
driven defect evolution, and processes in which diffusing 
solute atoms interact with mobile and/or immobile de- 
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feet structures. These, in practice, are inaccessible to 
conventional MD simulations. Conversely, PFC mod- 
els cannot reliably access the rapid time scales of MD 
and likely cannot predict some details of more complex 
defect atomic core structures as accurately as MD. Mi- 
croscopic phase field (MPF) models of defects^ are re- 
lated to PFC in the sense that both employ phase field 
methodologies, but MPF models utilize a top-down ap- 
proach, explicitly building the defect physics in by hand. 
PFC models employ a far simpler, atomic-level free en- 
ergy functional, from which all defect physics automat- 
ically emerge with fewer imposed constraints and with 
straightforward dynamical extensions. Nonetheless, by 
directly inco rporating ah initio data such as 7-surface 
energie d 49 * 50 !, MPF models have proven highly effective 
in terms of predicting static core-level defect structures 
and energies with greater quantitative accuracy than cur- 
rent PFC descriptions. At a more coarse-grained level, 
both discrete dislocation dynamics models 5 ^ and coarse- 
grained phase field dislocation dynamics models 48 read- 
ily describe length scales inaccessible to MD, MPF, and 
PFC. Similarly large length scales are in principle also 
accessible to coarse-grained complex amplitude represen- 



tations of PFC model^HHl^ which still retain atomistic 
resolution. 

Targeted PFC studies, examining more complex prob- 
lems in defect physics that involve atomic length scales 
and diffusive time scales, will be the subject of an upcom- 
ing publication. Topics that have been analyzed and will 
be addressed include climb fundamentals, the structure 
of jogged dislocation lines, jog constriction, jog pair anni- 
hilation, jog drag, screw dislocation cross-slip, formation 
and collapse of stacking fault tetrahedra from triangular 
Frank vacancy loops, Lomer-Cottrell or stair-rod disloca- 
tions, the interaction of dissociated edge and screw dis- 
locations with stacking fault tetrahedra, and dislocation 
creation mechanisms. 



Acknowledgments 

This work has been supported by the Natural Science 
and Engineering Research Council of Canada (NSERC), 
and access to supercomputing resources has been pro- 
vided by CLUMEQ/Compute Canada. 



1 K. R. Elder, M. Katakowski, M. Haataja, and M. Grant, 
Phys. Rev. Lett. 88, 245701 (2002). 

2 K. R. Elder and M. Grant, Phys. Rev. E 70, 051605 (2004). 

3 K. R. Elder, N. Provatas, J. Berry, P. Stefanovic, and 
M. Grant, Phys. Rev. B 75, 064107 (2007). 

4 N. Goldenfeld, B. P. Athreya, and J. A. Dantzig, Phys. 
Rev. E 72, 020601 (2005). 

5 B. P. Athreya, N. Goldenfeld, J. A. Dantzig, M. Green- 
wood, and N. Provatas, Phys. Rev. E 76, 056706 (2007). 

6 D. Yeon, Z. Huang, K. R. Elder, and K. Thornton, Philos. 
Mag. 90, 237 (2010). 

7 Z.-F. Huang, K. R. Elder, and N. Provatas, Phys. Rev. E 
82, 021605 (2010). 

8 R. Spatschek and A. Karma, Phys. Rev. B 81, 214201 
(2010). 

9 J. Berry, M. Grant, and K. R. Elder, Phys. Rev. E 73, 
031609 (2006). 

10 P. Stefanovic, M. Haataja, and N. Provatas, Phys. Rev. 
Lett. 96, 225504 (2006). 

11 P. Stefanovic, M. Haataja, and N. Provatas, Phys. Rev. E 

80, 046107 (2009). 

12 T. Hirouchi, T. Takaki, and Y. Tomita, Comp. Mat. Sci. 
44, 1192 (2009). 

13 P. Y. Chan, G. Tsekenis, J. Dantzig, K. A. Dahmen, and 
N. Goldenfeld, Phys. Rev. Lett. 105, 015502 (2010). 

14 A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, 
Phys. Rev. E 80, 031602 (2009). 

15 R. Backofen, A. Voigt, and T. Witkowski, Phys. Rev. E 

81, 025701(R) (2010). 

16 A. Jaatinen, C. V. Achim, K. R. Elder, and T. Ala-Nissila, 
Technische Mechanik 30, 169 (2010). 

17 J. Berry, K. R. Elder, and M. Grant, Phys. Rev. B 77, 
224114 (2008). 

18 J. Mellenthin, A. Karma, and M. Plapp, Phys. Rev. B 78, 



184110 (2008). 

19 D. L. Olmsted, D. Buta, A. Adland, S. M. Foiles, M. Asta, 
and A. Karma, Phys. Rev. Lett. 106, 046101 (2011). 

20 Y. Singh, Phys. Rep. 207, 351 (1991). 

21 D. Oxtoby, Annu. Rev. Mater. Res. 32, 39 (2002). 

22 T. V. Ramakrishnan and M. Yussouff, Phys. Rev. B 19, 
2775 (1979). 

23 The latter implies that the Fourier transformed kernel 
Ciik) is effectively shifted by a constant uuq relative to 
the kernel defined in Eq. |l]) . To maintain consistency with 
existing notation we will neglect this constant, and it will 
be implied that uhq must be subtracted from all kernels 

associated with Eq. ^ if one wishes to recover the 
kernel that is fully consistent with Eq. |l]). 

24 C. Dasgupta, Europhys. Lett. 20, 131 (1992). 

25 M. Greenwood, N. Provatas, and J. Rottler, Phys. Rev. 
Lett. 105, 045702 (2010). 

26 M. Greenwood, J. Rottler, and N. Provatas, Phys. Rev. E 
83, 031601 (2011). 

27 M. Greenwood, N. Ofori-Opoku, J. Rottler, and 
N. Provatas, Phys. Rev. B 84, 064104 (2011). 

28 J. K. Percus and G. J. Yevick, Phys. Rev. 110, 1 (1958). 

29 M. S. Wertheim, Phys. Rev. Lett. 10, 321 (1963). 

30 The {111} reflection will be broadened as well because the 
three additional {111} planes that are not parallel to the 
fault will contain l/(k z — ki^) 2 modes. 

31 Depending on local defect configurations some disloca- 
tion and grain boundary structures can annihilate by glide 
alone, leading to much lower energy barriers for removal. 
This effect is naturally captured by CDFT and PFC de- 
scriptions. 

32 S. A. Brazovskii, Sov. Phys.-JETP 41, 85 (1975). 

33 A. Jaatinen and T. Ala-Nissila, J. Phys.: Condens. Matter 
22, 205402 (2010). 



14 



G. I. Toth, G. Tegze, T. Pusztai, G. Toth, and L. Granasy, 
J. Phys.: Condens. Matter 22, 364101 (2010). 
K.-A. Wu, A. Adland, and A. Karma, Phys. Rev. E 81, 
061601 (2010). 

P. Y. Chan, N. Goldenfeld, and J. Dantzig, Phys. Rev. E 
79, 035701(R) (2009). 

J. Berry and M. Grant, Phys. Rev. Lett. 106, 175702 
(2011). 

J. P. Hirth and J. Lothe, Theory of Dislocations (John 
Wiley & Sons, Inc., 1982), 2nd ed. 

The independent elastic constants were measured numeri- 
cally by using Eq. (|3| and applying the appropriate strains, 
giving Cn = 0.0967, d 2 = 0.0388, and Caa — \i — 0.032. 
An isotropic approximation of Poisson's ratio is therefore 
v — Cnl{C\\ + C12) — 0.2863, and the anisotropy ra- 
tio is 2Caa/{Ch — C12) = 1.11. Thus there is a small de- 
gree of elastic anisotropy in this system, as evidenced by 
the slightly lower slopes for the anisotropic predictions in 
Fig. [8] The expressions for the Voigt and Reuss average 
anisotropic moduli can be found in Ref. 1381 
T. Byun, Acta Materialia 51, 3063 (2003). 



41 D. L. Olmsted, K. Y. Hardikar, and R. Phillips, Modelling 
and Simulation in Materials Science and Engineering 9, 
215 (2001). 

42 M. S. Duesbery, Modelling and Simulation in Materials 
Science and Engineering 6, 35 (1998). 

43 Z. S. Basinski, Phil. Mag. 4, 393 (1959). 

44 M. Bujard, G. Gremaud, and W. Benoit, J. Appl. Phys. 
62, 3173 (1987). 

45 F. R. N. Nabarro, Phil. Mag. A 75, 703 (1997). 

46 D. Mordehai, Y. Ashkenazy, I. Kelson, and G. Makov, 
Phys. Rev. B 67, 024112 (2003). 

47 D. Mordehai, I. Kelson, and G. Makov, Phys. Rev. B 74, 
184115 (2006). 

48 Y. Wang and J. Li, Acta Materialia 58, 1212 (2010). 

49 C. Shen and Y. Wang, Acta Materialia 52, 683 (2004). 

50 A. Hunter, I. J. Beyerlein, T. C. Germann, M. Koslowski, 
Phys. Rev. B 84, 144108 (2011). 

51 S. Groh and H. M. Zbib, Journal of Engineering Materials 
and Technology 131, 041209 (2009). 



